!> 来源：https://people.sc.fsu.edu/~jburkardt/f77_src/fftw/fftw.html
!> 本示例的许可证应该是：LGPL 许可证
program main

!*********************************************************************72
!
!  main is the main program for fftw_prb.
!
!  discussion:
!
!    fftw_prb tests the fftw library.
!
!  modified:
!
!    05 november 2007
!
!  author:
!
!    john burkardt
!
    implicit none

    call timestamp()
    write (*, '(a)') ' '
    write (*, '(a)') 'fftw_prb'
    write (*, '(a)') '  fortran77 version'
    write (*, '(a)') '  test the fftw library.'

    call test01()
    call test02()
    call test03()
    call test04()
!
!  terminate.
!
    write (*, '(a)') ' '
    write (*, '(a)') 'fftw_prb'
    write (*, '(a)') '  normal end of execution.'
    write (*, '(a)') ' '
    call timestamp()

end program main

subroutine test01()

!*********************************************************************72
!
!c test01: complex 1d data.
!
!  modified:
!
!    26 october 2005
!
    use fftw3_m
    implicit none

    integer n
    parameter(n=100)

    double precision a
    double precision b
    double precision r8_uniform_01
    integer i
    double complex in(n)
    double complex in2(n)
    double complex out(n)
    integer*8 plan_backward
    integer*8 plan_forward
    integer seed

    seed = 123456789

    write (*, '(a)') ' '
    write (*, '(a)') 'test01'
    write (*, '(a)') '  demonstrate fftw on a single vector '
    write (*, '(a)') '  of complex data.'
    write (*, '(a)') ' '
    write (*, '(a)') '  transform data to fft coefficients.'
    write (*, '(a)') '  backtransform fft coefficients to recover'
    write (*, '(a)') '  the data.'
    write (*, '(a)') '  compare recovered data to original data.'
!
!  compute the data.
!
    do i = 1, n
        a = r8_uniform_01(seed)
        b = r8_uniform_01(seed)
        in(i) = dcmplx(a, b)
    end do

    write (*, '(a)') ' '
    write (*, '(a)') '  input data:'
    write (*, '(a)') ' '

    do i = 1, n
        write (*, '(2x,i4,2x,2g14.6)') i, in(i)
    end do
!
!  make a plan for the fft, and forward transform the data.
!
    call dfftw_plan_dft_1d_(plan_forward, n, in, out, fftw_forward, fftw_estimate)

    call dfftw_execute_(plan_forward)

    write (*, '(a)') ' '
    write (*, '(a)') '  output fft coefficients:'
    write (*, '(a)') ' '

    do i = 1, n
        write (*, '(2x,i4,2x,2g14.6)') i, out(i)
    end do
!
!  make a plan for the backward fft, and recover the original data.
!
    call dfftw_plan_dft_1d_(plan_backward, n, out, in2, fftw_backward, fftw_estimate)

    call dfftw_execute_(plan_backward)

    write (*, '(a)') ' '
    write (*, '(a)') '  recovered input data divided by n:'
    write (*, '(a)') ' '

    do i = 1, n
        write (*, '(2x,i4,2x,2g14.6)') i, in2(i)/dble(n)
    end do
!
!  discard the information associated with the plans.
!
    call dfftw_destroy_plan_(plan_forward)
    call dfftw_destroy_plan_(plan_backward)

    return
end subroutine test01
subroutine test02()

!*********************************************************************72
!
!c test02: real 1d data.
!
!  modified:
!
!    23 october 2005
!
    use fftw3_m
    implicit none

    integer n
    integer nc

    parameter(n=100)
    parameter(nc=51)

    double precision r8_uniform_01
    integer i
    double precision in(n)
    double precision in2(n)
    double complex out(nc)
    integer*8 plan_backward
    integer*8 plan_forward
    integer seed

    seed = 123456789

    write (*, '(a)') ' '
    write (*, '(a)') 'test02'
    write (*, '(a)') '  demonstrate fftw on a single vector'
    write (*, '(a)') '  of real data.'
    write (*, '(a)') ' '
    write (*, '(a)') '  transform data to fft coefficients.'
    write (*, '(a)') '  backtransform fft coefficients to recover '
    write (*, '(a)') '  data.'
    write (*, '(a)') '  compare recovered data to original data.'
!
!  set up the input data, a real vector of length n.
!
    do i = 1, n
        in(i) = r8_uniform_01(seed)
    end do

    write (*, '(a)') ' '
    write (*, '(a)') '  input data:'
    write (*, '(a)') ' '

    do i = 1, n
        write (*, '(2x,i4,2x,g14.6)') i, in(i)
    end do
!
!  set up a plan, and execute the plan to transform the in data to
!  the out fft coefficients.
!
    call dfftw_plan_dft_r2c_1d_(plan_forward, n, in, out, fftw_estimate)

    call dfftw_execute_(plan_forward)

    write (*, '(a)') ' '
    write (*, '(a)') '  output fft coefficients:'
    write (*, '(a)') ' '

    do i = 1, nc
        write (*, '(2x,i4,2x,2g14.6)') i, out(i)
    end do
!
!  set up a plan, and execute the plan to backtransform the
!  complex fft coefficients in out to real data.
!
    call dfftw_plan_dft_c2r_1d_(plan_backward, n, out, in2, fftw_estimate)

    call dfftw_execute_(plan_backward)

    write (*, '(a)') ' '
    write (*, '(a)') '  recovered input data divide by n:'
    write (*, '(a)') ' '

    do i = 1, n
        write (*, '(2x,i4,2x,g14.6)') i, in2(i)/dble(n)
    end do
!
!  release the memory associated with the plans.
!
    call dfftw_destroy_plan_(plan_forward)
    call dfftw_destroy_plan_(plan_backward)

    return
end subroutine test02
subroutine test03()

!*********************************************************************72
!
!c test03: complex 2d data.
!
!  modified:
!
!    05 november 2007
!
!  author:
!
!    john burkardt
!
    use fftw3_m
    implicit none

    integer nx
    parameter(nx=8)
    integer ny
    parameter(ny=10)

    double precision a
    double precision b
    double precision r8_uniform_01
    integer i
    double complex in(nx, ny)
    double complex in2(nx, ny)
    integer j
    double complex out(nx, ny)
    integer*8 plan_backward
    integer*8 plan_forward
    integer seed

    seed = 123456789

    write (*, '(a)') ' '
    write (*, '(a)') 'test03'
    write (*, '(a)') '  demonstrate fftw on a 2d complex array'
    write (*, '(a,i8)') '  nx = ', nx
    write (*, '(a,i8)') '  ny = ', ny
    write (*, '(a)') ' '
    write (*, '(a)') '  transform data to fft coefficients.'
    write (*, '(a)') '  backtransform fft coefficients to recover'
    write (*, '(a)') '  the data.'
    write (*, '(a)') '  compare recovered data to original data.'
!
!  compute the data.
!
    do j = 1, ny
        do i = 1, nx
            a = r8_uniform_01(seed)
            b = r8_uniform_01(seed)
            in(i, j) = dcmplx(a, b)
        end do
    end do

    write (*, '(a)') ' '
    write (*, '(a)') '  input data:'
    write (*, '(a)') ' '

    do i = 1, nx
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,2g14.6)') i, j, in(i, j)
        end do
    end do
!
!  make a plan for the fft, and forward transform the data.
!
    call dfftw_plan_dft_2d_(plan_forward, nx, ny, in, out, fftw_forward, fftw_estimate)

    call dfftw_execute_(plan_forward)

    write (*, '(a)') ' '
    write (*, '(a)') '  output fft coefficients:'
    write (*, '(a)') ' '

    do i = 1, nx
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,2g14.6)') i, j, out(i, j)
        end do
    end do
!
!  make a plan for the backward fft, and recover the original data.
!
    call dfftw_plan_dft_2d_(plan_backward, nx, ny, out, in2, fftw_backward, fftw_estimate)

    call dfftw_execute_(plan_backward)

    write (*, '(a)') ' '
    write (*, '(a)') '  recovered input data divided by nx * ny:'
    write (*, '(a)') ' '

    do i = 1, nx
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,2g14.6)') i, j, in2(i, j)/dble(nx*ny)
        end do
    end do
!
!  discard the information associated with the plans.
!
    call dfftw_destroy_plan_(plan_forward)
    call dfftw_destroy_plan_(plan_backward)

    return
end subroutine test03
subroutine test04()

!*********************************************************************72
!
!c test04: real 2d data.
!
!  discussion:
!
!    in contrast to the c example, in fortran it is the first dimension
!    of the complex coefficient array that is "half" the size.
!
!  modified:
!
!    05 november 2007
!
!  author:
!
!    john burkardt
!
    use fftw3_m
    implicit none

    integer nx
    parameter(nx=8)
    integer ny
    parameter(ny=10)

    integer nh
    parameter(nh=(nx/2) + 1)

    double precision a
    double precision b
    double precision r8_uniform_01
    integer i
    double precision in(nx, ny)
    double precision in2(nx, ny)
    integer j
    double complex out(nh, ny)
    integer*8 plan_backward
    integer*8 plan_forward
    integer seed

    seed = 123456789

    write (*, '(a)') ' '
    write (*, '(a)') 'test04'
    write (*, '(a)') '  demonstrate fftw on a 2d real array'
    write (*, '(a,i8)') '  nx = ', nx
    write (*, '(a,i8)') '  ny = ', ny
    write (*, '(a)') ' '
    write (*, '(a)') '  transform data to fft coefficients.'
    write (*, '(a)') '  backtransform fft coefficients to recover'
    write (*, '(a)') '  the data.'
    write (*, '(a)') '  compare recovered data to original data.'
!
!  compute the data.
!
    do j = 1, ny
        do i = 1, nx
            in(i, j) = r8_uniform_01(seed)
        end do
    end do

    write (*, '(a)') ' '
    write (*, '(a)') '  input data:'
    write (*, '(a)') ' '

    do i = 1, nx
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,g14.6)') i, j, in(i, j)
        end do
    end do
!
!  make a plan for the fft, and forward transform the data.
!
    call dfftw_plan_dft_r2c_2d_(plan_forward, nx, ny, in, out, fftw_estimate)

    call dfftw_execute_(plan_forward)

    write (*, '(a)') ' '
    write (*, '(a)') '  output fft coefficients:'
    write (*, '(a)') ' '

    do i = 1, nh
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,2g14.6)') i, j, out(i, j)
        end do
    end do
!
!  make a plan for the backward fft, and recover the original data.
!
    call dfftw_plan_dft_c2r_2d_(plan_backward, nx, ny, out, in2, fftw_estimate)

    call dfftw_execute_(plan_backward)

    write (*, '(a)') ' '
    write (*, '(a)') '  recovered input data divided by nx * ny:'
    write (*, '(a)') ' '

    do i = 1, nx
        do j = 1, ny
            write (*, '(2x,i4,2x,i4,2x,g14.6)') i, j, in2(i, j)/dble(nx*ny)
        end do
    end do
!
!  discard the information associated with the plans.
!
    call dfftw_destroy_plan_(plan_forward)
    call dfftw_destroy_plan_(plan_backward)

    return
end subroutine test04
function r8_uniform_01(seed)

!*********************************************************************72
!
!c r8_uniform_01 returns a unit double precision pseudorandom number.
!
!  discussion:
!
!    this routine implements the recursion
!
!      seed = 16807 * seed mod ( 2**31 - 1 )
!      r8_uniform_01 = seed / ( 2**31 - 1 )
!
!    the integer arithmetic never requires more than 32 bits,
!    including a sign bit.
!
!    if the initial seed is 12345, then the first three computations are
!
!      input     output      r8_uniform_01
!      seed      seed
!
!         12345   207482415  0.096616
!     207482415  1790989824  0.833995
!    1790989824  2035175616  0.947702
!
!  modified:
!
!    23 october 2005
!
!  author:
!
!    john burkardt
!
!  reference:
!
!    paul bratley, bennett fox, l e schrage,
!    a guide to simulation,
!    springer verlag, pages 201-202, 1983.
!
!    pierre l'ecuyer,
!    random number generation,
!    in handbook of simulation,
!    edited by jerry banks,
!    wiley interscience, page 95, 1998.
!
!    bennett fox,
!    algorithm 647:
!    implementation and relative efficiency of quasirandom
!    sequence generators,
!    acm transactions on mathematical software,
!    volume 12, number 4, pages 362-376, 1986.
!
!    p a lewis, a s goodman, j m miller,
!    a pseudo-random number generator for the system/360,
!    ibm systems journal,
!    volume 8, pages 136-143, 1969.
!
!  parameters:
!
!    input/output, integer seed, the "seed" value, which should not be 0.
!    on output, seed has been updated.
!
!    output, double precision r8_uniform_01, a new pseudorandom variate,
!    strictly between 0 and 1.
!
    use fftw3_m
    implicit none

    double precision r8_uniform_01
    integer k
    integer seed

    k = seed/127773

    seed = 16807*(seed - k*127773) - k*2836

    if (seed < 0) then
        seed = seed + 2147483647
    end if
!
!  although seed can be represented exactly as a 32 bit integer,
!  it generally cannot be represented exactly as a 32 bit real number!
!
    r8_uniform_01 = dble(seed)*4.656612875d-10

    return
end function r8_uniform_01
subroutine timestamp()

!*********************************************************************72
!
!c timestamp prints out the current ymdhms date as a timestamp.
!
!  discussion:
!
!    this fortran77 version is made available for cases where the
!    fortran90 version cannot be used.
!
!  modified:
!
!    16 september 2005
!
!  author:
!
!    john burkardt
!
!  parameters:
!
!    none
!
    use fftw3_m
    implicit none

    character(len=8) date
    character(len=10) time

    call date_and_time(date, time)

    write (*, '(a8,2x,a10)') date, time

    return
end subroutine timestamp
